GENXICC2.1: An Improved Version of GENXICC for Hadronic Production of 

Doubly Heavy Baryons 
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We present an improved version of GENXICC, which is a generator for hadronic production of 
the doubly heavy baryons H cc , Sf, c and S&f, and has been raised by C.H. Chang, J.X. Wang and 
X.G Wu [Comput. Phys. Commun. 177 (2007) 467; Comput. Phys. Commun. 181 (2010) 
1144]. In comparison with the previous GENXICC versions, we update the program in order to 
generate the unweighted baryon events more effectively under various simulation environments, 
whose distributions are now generated according to the probability proportional to the integrand. 
One Les Houches Event (LHE) common block has been added to produce a standard LHE data 
file that contains useful information of the doubly heavy baryon and its accompanying partons. 
Such LHE data can be conveniently imported into PYTHIA to do further hadronization and decay 
simulation, especially, the color-flow problem can be solved with PYTHIA8.0. 
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NEW VERSION PROGRAM SUMMARY 

Title of program : GENXICC2.1 

Program obtained from : CPC Program Library 

Reference to original program : GENXICC 

Reference in CPC : Comput. Phys. Commun. 177, 467 
(2007); Comput. Phys. Commun. 181, 1144 (2010) 

Does the new version supersede the old program: No 

Computer : Any LINUX based on PC with FORTRAN 
77 or FORTRAN 90 and GNU C compiler as well 

Operating systems : LINUX 

Programming language used : FORTRAN 77/90 

Memory required to execute with typical data : About 
2.0 MB 

No. of bytes in distributed program : About 2 MB, 
including PYTHIA6.4 

Distribution format : .tar.gz 

Nature of physical problem : Hadronic production of 
doubly heavy baryons S cc , S;, c and 'Ebb- 

Method of solution : The upgraded version with proper 
interface to PYTHIA can generate full production 



and decay events, either weighted or unweighted, con- 
veniently and effectively. Especially, the unweighted 
events are generated by using an improved hit-and-miss 
approach. 

Reasons for new version : Responding to the feed- 
back from users of CMS and LHCb groups at the large 
hadronic collider, and basing on the recent improvements 
of PYTHIA on the color-flow problem, we improve the 
efficiency for generating the unweighted events, and also 
improve the color-flow part for further hadronization. 
Especially, an interface has been added to import 
the output production events into a suitable form for 
PYTHIA8.0 simulation, in which the color-flow during 
the simulation can be correctly set. 

Typical running time : It depends on which option is 
chosen to match PYTHIA when generating the full 
events and also on which mechanism is chosen to gener- 
ate the events. Typically, for the dominant gluon-gluon 
fusion mechanism to generate the mixed events via 
the intermediate diquarks in (cc)[ 3 5i]3 and (cc)[ 1 S'o]6 
states, setting IDWTUP=3 and unwght=.true., it takes 
30 minutes to generate 10 5 unweighted events on a 
2.27GHz Intel Xeon E5520 processor machine; setting 
IDWTUP=3 and unwght=. false, or IDWTUP=1 and 
IGENERATE=0, it only needs 2 minutes to generate 
the 10 5 baryon events (the fastest way, for theoretical 
purpose only). As a comparison, for previous GENXICC 
versions, if setting IDWTUP^l and IGENERATE=1, it 
takes about 22 hours to generate 1000 unweighted events. 

Keywords : Event generator; Doubly heavy baryons; 
Hadronic production. 
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Summary of the changes (improvements) : 1) The scheme 
for generating unweighted events has been improved; 2) 
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One Les Houches Event (LHE) common block has been 
added to record the standard LHE data in order to be 
the correct input for PYTHIA8.0 for later simulation; 
3) We present the code for connecting GENXICC to 
PYTHIA8.0, where three color-flows have to be correctly 
set for later simulation. More specifically, we present the 
changes together with their detailed explanations in the 
following : 

• Unweighted events generation. For theoreti- 
cal studies, e.g. to derive the total baryon pro- 
duction cross-section or various differential distri- 
butions, one can directly use the fastest way, e.g. 
setting the PYTHIA parameter IDWTUP=3 and 
unwght=.false. or setting IDWTUP= 1 and IGEN- 
ERATE=0 (in these cases, xmaxup should be set 
as 0), to generate the baryon events [l|. By using 
GENXICC @, H[ in this way, some interesting prop- 
erties for hadronic production of 5 CC , S(, c and Sfcj, 
have been found in the literature, cf. Refs.jjHq]. 
While, for the events simulation in detector condi- 
tions, it is necessary to get the unweighted events. 
In previous GENXICC versions, the unweighted 
events are generated by setting IDWGTUP = 1 
and IGENERATE= 1; i.e., the events are gen- 
erated according to PYTHIA's inner mechanism, 
the so-called hit-and-miss approach (von Neumann 
algorithm), to reject those unsatisfied events and 
output the allowed events. But, as is well-known, 
the original hit-and-miss approach is really time- 
consuming. Some alterations must be made to im- 
prove its efficiency. 

As an intermediate step, in BCVEGPY2.1a Q 
we have suggested a practical trick to increase 
the efficiency of generating unweighted events 
(BCVEGPY is a generator for hadronic production 
B c meson Q). In this trick, other than choosing 
the maximum differential cross-section as a refer- 
ence weight in the hit-and-miss approach, we di- 
rectly select an effective differential cross-section, 
which is smaller than the maximum one, as the 
reference weight Q. This treatment can largely 
improve the generation efficiency without affecting 
the total cross-section of the process. However, in 
using this trick to generate unweighted events such 
as for CMS detector simulation, one will inciden- 
tally find a false peak in the B c -p t distributions. 
This is caused by the fact that sometimes the same 
event will be stored with (false) large number of 
times in the hit-and-miss process. Then, we are 



facing a dilemma: such a false peak can be avoided 
by rising the effective reference weight to a value 
approaching the maximum weight, but, inversely, a 
larger reference weight will surely lead to a much 
longer running-time. 

One observes that by using the VEGAS algo- 
rithm @, the SPRING-BASES program per- 
forms the integration in using the BASES subrou- 
tines and generates events with a probability pro- 
portional to the integrand in using the SPRING 
subroutines. After each iteration of VEGAS run- 
ning, the integration result and the maximum value 
of the function will be stored in a file for each cell 
of the adaptive mesh. In the generation stage, a 
cell is chosen with a probability proportional to 
the corresponding integral, and then a point in the 
cell is generated using the hit-and-miss approach. 
This method is highly efficient, but it has the disad- 
vantage that the required amount of storage space 
grows exponentially with the integration dimen- 
sion. 

Next, in POWHEG program [ll[ the authors have 
developed a new method MINT [12j to replace 
the SPRING-BASES package. This MINT package 
also use the VEGAS algorithm to perform the in- 
tegration. What's the difference is that it does not 
store the value of the integral but stores the upper 
bound value for each cell. The multidimensional 
stepwise function that equals to the upper bound 
of the function to be integrated in each cell is in 
fact an upper bound for the whole function, which 
is the wanted upper bound for BCVEGPY2.1a or 
the PYTHIA. So, the program is to find the up- 
per bound grid for those cells. And next, by using 
again the hit-and-miss technique in each cell, one 
can generate the points according to the original 
distribution. 

Basing on these methods, as a further improve- 
ment, we present an ultimate solution to gener- 
ate unweighted events in the present new GENX- 
ICC version. We adopt the MINT algorithm but 
with certain alterations to do the simulation. For 
the purpose, we change the VEGAS subroutine as 
follows. Three new variables have been added in 
the original VEGAS subroutine, where xint is the 
integral value for the integrand f xn after a ndim- 
dimensional integration, the xmax array records the 
upper bounding envelope of the integrand in all 
cells, imode is a flag : 



vegas (f xn , ndim , ncall , itmx , nprn , xint , xmax , imode ) 



When called with imode=0, vegas performs 
the integration over the integrand fxn, and 
stores the answer in a common block parame- 



ter vegsec. 

— xmax stands for a (nvegbin,ndim) dimensional 
array, where nvegbin denotes the bin num- 
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ber for each coordinate, ndim stand for the 
integration dimension. When called with 
imode=l, vegas will first initiate all the ele- 
ments of xmax to be xint 1/ndim , where xint 
equals to the value of vegsec that has been 
derived from a previous VEGAS running with 
imode=0. During the following sampling it- 
eration, when the calculated integral value 
is larger than the initial xmax(nvegbin,ndim) 
value in a specific cell, then the value of 
xmax(nvegbin,ndim) for this cell will be in- 
creased by a fixed factor / = 1 + 1/ 10 ndim. 
After a sufficiently large number of calls, the 
values of xmax(nvegbin,ndim) will be stabi- 
lized for all cells. Such a final xmax array will 
be stored in the same grid file as that of the 
importance sampling function in order to do 
the final simulation. 
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FIG. 1. Comparison of the normalized S cc transverse mo- 
mentum (Pt) and rapidity (y) distributions derived by 
IDWTUP=3 (events) and IDWTUP=1 (differential cross- 
sections), which are represented by solid line and dotted line 
respectively. 



Comparing to the previous GENXICC versions, in 
doing the initialization (subroutine evntinit), we 
will call vegas twice with imode=0 and imode=l 
accordingly to generate the upper bound grid xmax 
and also a more precise importance sampling func- 
tion. Practically, the user can directly use the 
existed grid file derived by previous VEGAS run- 
ning to generate events by setting methodevnt=2 
or methodevnt=3 without running VEGAS again, 
which is the same as the older GENXICC versions. 

Once the xmax array has been set up in previous 
steps, one can call the subroutine gen to generate 
events. For the purpose, three options for calling 
gen subroutine are programmed : 

jmode=0 

call gen(fxn, ndim, xmax, jmode) 

jmode=l 

do j=l, 10000 

call gen(fxn, ndim, xmax, jmode) 

enddo 
jmode=3 

call gen(fxn, ndim, xmax, jmode) 

where jmode=0 is to initializes a step- wise function 
xmmm which descripted infill]. And jmode=3 is to 
print out the generation statistics. 

The calling for the gen subroutine with jmode=l 
has been implemented into the UPEVNT subrou- 
tine to generate events according to the probabil- 
ity proportional to the integrand. Each event pro- 
duced needs several times of iteration with three 
steps procedure as follows: 

1. Calculate upper bounding function by gener- 
ating a set of step- wise functions, each of them 
is associated with a specific coordinate (di- 
mension) . 



2. Call the phase_gen subroutine to generate a 
random phase-space point and calculate the 
integral. 

3. Judge whether such point be kept or not by 
using the hit-and-miss approach with the help 
of the upper bounding function. 

In VEGAS the integral together with its numerical 
error are related to the sampling numbers ncall 
and the iteration times itmx. So, to generate full 
events, we suggest the user to do a test running 
first in order to find an effective and time-saving 
parameters for VEGAS. Furthermore, to validate 
the program, we use the same default parameters 
as the input for the program to generate mixed 
events via the intermediate diquark in (cc)[ 3 Si}^ 
and (cc)[ 1 S'o]6 states, and the same for other two 
doubly heavy baryons Sf, c and S^f,. 

As a cross-check, we derive the unweighted S cc 
event distributions by setting IDWTUP=3 and 
unwght=.true., and the weighted S cc differential 
distributions by setting IDWTUP=1 and IGEN- 
ERATE=0, respectively, which are shown in FIGQ] 
The two distributions after proper normalization 
agree well with each other. This demonstrates that 
our present scheme for unweighted events is correct. 

• Color-flow problem. Within the framework of 
non-relativistic QCD (NRQCD) [lj], the produc- 
tion of Sqq/ baryon can be factorized into two 
steps: The first step is to produce two free heavy- 
quark pairs QQ and Q'Q' , which is perturbatively 
calculable. The second step is to make the two 
heavy quarks Q and Q' into a bounding diquark 
(QQ') in [ 3 5i] (or [ l S ]) spin state and in 3 (or 6) 
color state accordingly; then it will be hadronized 
into Sqq< baryon by grabbing a light quark u or d 
or s (plus suitable number of gluons) , whose prob- 
ability is described by non-perturbative NRQCD 
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matrix element. More explicitly, the intermedi- 
ate diquarks in S cc and S&6 have two spin-and- 
color configurations [ 3 Si]3 and ["""Sole! while for 
the intermediate diquark (be) in Eft c , there are four 
spin-and-color conhgurations Sb c [ 3 Si]g, Sb c [ 3 Si]6, 
^■bc^So}^, and Eb^ 1 Sq]s- 

Since a baryon is constructed by three valance 
quarks, under the standard color-flow decomposi- 
tion, there must be three different color-flow lines 
being ended at a baryon 0, [l4| ■ It is different from 
the case of meson, where the color-flow lines of the 
quark and anti-quark inside a meson are contin- 
ued. The previous PYTHIA6.4 can only generate 
full events with two or less independent color-flow 
lines, thus in GENXICC2.0, we adopt a 'cheating 
method' to generate the events. That is, by using 
the fact that 3(g)3 = 6®3and3(g)3 = 801in 
general QCD SU(3) color space [3j : 

— We combine any two of the color-flow lines 
ended with two quarks into one anti-color-flow 
line ended with one anti-quark with a color 3 
that is different from the two quarks (the third 
color in respect to those of the two quarks); 

— Secondly, such anti-color-flow line obtained 
by the combination may be continued (con- 
nected) to the remaining quark's color-flow 
line in the baryon; 

— Finally, as a consequence, the color-flow lines 
ended at a baryon become 'joined without 
ends' at all, which is the requirement of the 
color-singlet bound state. 

However we should point out that due to approxi- 
mation and simplification with 'cheating method', 
the obtained information about the 'tiny jets', cor- 
responding to the soft anti-quark and soft gluon(s) 
produced in fragmentation of doubly heavy di- 
quark, may not be very reliable. When the ex- 
periment analyzer uses the generator to simulate 
the baryon decay and other parton hardronization, 
they are still facing the color-flow rearrangement 
error; sometimes, PYTHIA will present an error 
message to show that the color-flow rearrangement 
is wrong during the parton's evolution process, and 
then it will stop running. 

To generate full events of the doubly heavy baryons 
smoothly, the best way is to improve PYTHIA with 
proper treatment on the color-flow lines ended at 
the baryon. Fortunately, such an improvement has 



been done in its newest version PYTHIA8.0. Based 
on the suggestion from Peter Skands, we find that 
the further event simulation can be implemented 
into PYTHIA8.0 correctly even with the prev ious 
generated Les Houches Event (LHE) files [l5[. As 
has been described in Ref. [l6[ , the read-in of exter- 
nal generator's LHE files generated by PYTHIA6.4 
is simply technically less sophisticated and less able 
to deal with junctions, even though the physics im- 
plementation of junction fragmentation is in princi- 
ple the same. As an solution, PYTHIA8.0 improves 
the treatment on these LHE files. 

We adopt the same trick as that of 
BCVEGPY2.1a Q to generate and record 
the data; i.e. two subroutines have been in- 
troduced in the file pythialheinit.F. One of the 
subroutine XICC.PYUPIN is used to fill the 
HEPRUP common block with information on the 
incoming beams and the allowed processes, and 
optionally stores that information on file. Another 
subroutine XICC_WRITE_LHE is used to store 
event information in the HEPEUP common block. 
And these two subroutines are called by the main 
program xicc.F to generate the LHE file that 
records the momentum and color information for 
the events }17| . 

More specifically, in the main program, the sub- 
routine UPEVNT will be called for generating the 
baryon events, which is used to call the program 
to generate the baryon with a probability propor- 
tional to the importance sampling function. Here, 
one can also use the PYTHIA subroutine PYEVNT 
for the purpose, but one should at the same time 
switch off the hadronization, the initial and final 
state parton shower and so on, in order to avoid the 
color-flow rearrangement error. Those generated 
baryon information together with the information 
of the accompanying partons will be stored in the 
Les Houches common block and will be export to 
a LHE file "GENXICC.lhe". And then, such LHE 
file can be used when necessary by PYTHIA8.0 to 
do the following simulations. Here, to successfully 
simulate the baryon's production and decay, the 
user need to install the PYTHIA8[ll| following the 
instruction of its official web-site : 



http://home.thep.lu.se/~torbjorn/Pythia.html 



For using our generator, the user can use the fol- 
lowing command to compile the configuration file, 



g++ -02 -ansi -pedantic -W -Wall -Wshadow -1$ (PYTHIA8) /include genxicc.ee -o bin/genxicc . exe 
-L$(PYTHIA8) /lib/archive -lpythia8 -llhapdf dummy 

I 

where the $ (PYTHIA8) stands for the PYTHIA8.0 installation directory. 
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For convenience, we put an example configuration 
file in the package for generating the full the baryon 
production and decay events in PYTHIA8.0, which 
is placed in the main folder of the program and is 
named as "genxicc.cc" . 
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